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We propose a new recursive procedure to estimate the microcanonical density of states in mul- 
ticanonical Monte Carlo simulations which relies only on measurements of moments of the energy 
distribution, avoiding entirely the need for energy histograms. This method yields directly a piece- 
wise analytical approximation to the microcanonical inverse temperature, (3{E), and allows improved 
control over the statistics and efficiency of the simulations. We demonstrate its utility in connec- 
tion with recently proposed schemes for improving the efficiency of multicanonical sampling, either 
with adjustment of the asymptotic energy distribution or with the replacement of single spin flip 
dynamics with collective updates. 



I. INTRODUCTION 

Over the past decade, the use of Monte Carlo methods 
has broken the boundaries of statistical physics and 
has become a widely used computational tool in fields 
as diverse as chemistry, biology and even sociology or 
finance. 

Despite the enormous success of the well known 
Metropolis importance sampling algorithm, its narrow 
exploration of the phase space and characteristic con- 
vergence difficulties motivated a number of different ap- 
proaches. Firstly, the techniques for harvesting useful 
information from the statistical data obtained at a given 
temperature were improved in order to extrapolate the 
results in a small temperature range and, therefore, re- 
duce the number of required independent simulations 
[. Secondly, cluster update algorithms were proposed 
, 15 in order to overcome critical slowing down. Fi- 
nally, the requirement of constant temperature was lifted, 
allowing the system to explore a wider range of the en- 
ergy spectrum. In Simulated Tempering, the tempera- 
ture becomes a dynamical variable, which can change in 
the Markov process Parallel Tempering uses several 
replicas of the system running at different temperatures 
and introduces the swapping of configurations between 
the various Markov chains [a 0, E| • The flat histogram 
approach 0,^3 replaces the asymptotic Boltzmann dis- 
tribution for the asymptotic probability of sampling a 
given state of energy Ei, pi oc ex.'p[j3Ei), by pi (x l/n(E), 
n{E) being the microcanonical density of states, thus en- 
suring that every energy is sampled with equal probabil- 
ity. 

The first major obstacle to this last approach is the 



obvious difficulty of accessing the true density of states 
of a given system; several clever algorithms have been 
proposed to this end Nevertheless, even in cases 

where the true density of states is known a priori, recent 
studies 0, 01 have shown that very long equilibration 
times can remain a serious concern with multicanonical 
methods. Furthermore, the number of independent sam- 
ples is strongly dependent on energy, making error es- 
timation rather tricky T^. This has led some authors 
(El) to question the reliability of the results obtained 
by the Multicanonical Method in spin glass models at 
low temperatures. To address these issues, the require- 
ment of a perfectly flat histogram was also lifted [l8l |. 
sacrificing the equal probability of sampling each energy 
in favor of minimizing tunneling times. 

In section II, we propose a variation of the algorithm 
to estimate the density of states that does not use his- 
tograms, instead relying entirely of measurements of cu- 
mulants of the energy distribution at each stage of the 
simulation to build a piecewise analytical approximation 
to the statistical entropy, S{E) = \n{n{E)), and inverse 
temperature P{E) = dS{E)/dE. We find that the time 
required to explore the entire energy spectrum scales 
more favorably with system size than histogram based 
methods. The method is as easy to apply in systems 
with continuous spectrum as in the discrete case, can be 
quite naturally adapted to running a simulation in a cho- 
sen energy range, and accommodates without difficulty 
tunneling times optimization schemes and cluster update 
methods. 

In Section III we show that the optimization scheme 
proposed in can be applied during the process of 
estimating the density of states, still avoiding histograms, 
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and with significant efficiency gains. 

In Section IV, we demonstrate tlie usefulness of tlic 
analytical approximation of f3{E) in generalizing Wolff's 
cluster algorithm "B"] to multicanonical simulations. This 
generalization maintains an acceptance probability still 
very close to unity, growing large clusters at low temper- 
atures and small ones at higher temperatures. The opti- 
mization procedure reported in 18] is also implemented 
for this cluster dynamics. 

This algorithm has already been applied with suc- 
cess to both discrete (Ising models on regular lattices 
and Small World networks, Ising Spin Glasses) and 
continuous models (XY and Heisenberg models, with 
both short and long-range interactions, namely dipolar 
interactions) , but this is its first systematic presenta- 
tion. 



II. PROPOSAL 

The usual way of ensuring an asymptotic distribution 
Pi cx exp[—uj{Ei)] is to use a Markov chain algorithm 
in which the transition probability to go from state i to 
state j with respective energies Ei and Ej is 



Wij ~ min(l. 



-Aw(£;) 



(1) 



This choice leads to an asymptotic energy distribution 
probability given by 



E{E) cx exp(5(i;) (£;)), 



(2) 



where the entropy is defined as the logarithm of the (un- 
known) density of energy states, S(E) — ln(n(£')). In 
principle, this relation allows a calculation of the entropy 
S{E^ (up to an irrelevant constant So) from a numerical 
determination of the distribution H[E) for any value of 
the energy E simply as: 



S{E) ^Sq+lu{E)+ \tlH{E) 



(3) 



However, this is not always efficient for any choice of the 
function lj{E). Consider, for example, the widely used 
Metropolis choice, u{E) = l3E (with /? = l/T, the in- 
verse temperature). In many systems, the resulting dis- 
tribution H{E) has the shape of a bell curve, usually ap- 
proximated by a Gaussian distribution (Figure^. Since 
it is very unlikely to generate statistically significant con- 
figurations in the tails of the distribution for H{E), the 
usefulness of the formula ^ is limited to values of E not 
too far from the mean value /i/5. "Not too far" means ex- 
plicitly that the above formula is limited to those values 
of E such that \E — fi/^l < aap with a'^2 — 4:. The mean 
value /i/3 and the variance ct^ of the distribution H{E) 
satisfy: 



dS_ 
dE 



/9; 
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dE"^ 
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(4) 



A clever choice for uj{E) can greatly improve the range 
of values of E for which the formula © is useful. For 
instance, if we were to chose lo{E) cx hi{n{E)) then the 
resulting distribution H{E) would be constant in E and 
all the energy values would be sampled with the same 
frequency. However, it is clear that this is impossible 
since n(E) is precisely the function wc want to determine. 

Berg's scheme uses a series of functions u!i{E), each 
one of them gives information on n{E) for a range of en- 
ergy values. The initial choice ijJo{E) = (equivalent to 
a Metropolis choice at infinite temperature) provides a 
histogram Hq(E) from which one derives an estimation 
of the entropy, So{E), valid for those values of E visited 
in a statistically significant way. After this stage, a new 
simulation is performed with uJi{E) ~ Sq{E) in the re- 
gion visited in the previous simulation, from which we 
obtain an entropy estimate Si{E) valid in another range 
of energies, and so on. Berg proposes a recursion scheme 
that allows systematic corrections of uji (E) at all visited 
energies, and ensures the convergence to the true entropy 
for any energy, so ensuring a flat energy histogram. 
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Figure 1; Schematic representation of the true derivative of 
the entropy S{E) and two different choices for the deriva- 
tive of uj(E), with the corresponding energy probability dis- 
tributions. In panel a), the Metropolis approach, where 
ui{E) = /3E and, in panel b), the ideal extended ensemble 
approach with u}{E) = S{E). 

This procedure has some limitations. To begin with, 
the entropy can only be estimated inside the energy range 
visited in the last simulation and a bad choice for the 
entropy outside this region can severely limit the explo- 
ration of lower energies; secondly, rarely visited energies 
introduce a large error in the estimated entropy (hence 
the need for recursion introduced by Berg in order to min- 
imize this error). Furthermore, the fact that one counts 
visits in each energy implies that, for continuous systems, 
the energy spectrum must be discretized. 

Thermodynamic functions such as the entropy can be 
treated as continuous functions of energy, both for dis- 
crete and continuous spectra, for not too small systems. 
We make use of the Gaussian approximation: 

S{E)^S{fif,) + P{E-fif,)-^{E-fipf (5) 
to propose the sequence of weight functions uj(E). 
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Figure 2: In the current proposal, the weight function ijj{E) 
approximates the entropy between two energies and and 
is linear in energy outside this region, duj/dE = Pr for E < /ir 
and du/dE — f3o for E > /iq- 
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Figure 3: A comparison of the evolution of the multicanonical 
recursion in the present scheme (upper panel) and the one 
described in 21] (lower panel). The simulation was done for 
a Ising square sample with L = 64, and, for simplicity, we 
only show histograms from 1 out of each 5 runs. 



A. Comparison with histogram based recursion 



After running an initial simulation at a inverse tem- 
perature /3o, with luq = PqE, and measuring and ctq, 
we modify the weight function to lui(E) defined by: 

, / p^ f Po E > fj.0 

f3,-^{E-^io) m<i?</^o (6) 

with /ii = /io — ctcTo and Pi — Po + ci/o'o- In a new 
simulation with transition rates Wij = min(l, e"'^'^^'^^-') 
we obtain H{E) ~ constant, for energies such that /ii < 
E < Ho, while for E < fii, H{E) is a "half Gaussian" 
with maximum at fii. It is straightforward to estimate 
(Ti from the average of {E — fii) for E < ^i. We are 
therefore be able to add another branch to lo{E), and 
so on, until we reach the lowest temperature we wish to 
study. On the iteration of order r the histogram is flat 
between /io and and half Gaussian below /ir (Figure!^. 
Using u){E) = Pq for E > fiQ, we effectively restrict the 
simulation to energies below /ip, apart from a Gaussian 
tail above this energy. 

The fact that P{E) — const in the unexplored energy 
regions, corresponding to Boltzmann sampling, means 
that we can use all techniques developed for the Metropo- 
lis algorithm in order to be confident on the results ob- 
tained in this region, before moving on to lower energies. 
These regions of canonical sampling can also be used to 
restrict the simulation to a specific temperature range, 
for instance around a critical temperature, which can 
dramatically increase the efficiency and precision of the 
simulation. The method can also be refined by keeping 
higher order terms in the expansion of the entropy which 
can be obtained by measuring higher order moments of 
the energy. We have successfully used an expansion of 
S{E) up to fourth order terms. 



We compared an implementation of this proposal 
based on classical fiuctuation theory (CFP) with the mul- 
ticanonical recursion as described in _21| in a simulation 
of a two dimensional Ising model: 

H = — J SiSj. 

<ij> 

Except for the method of estimating the density of states 
(or rather, its logarithm, the entropy), the results were 
obtained using exactly the same code, specifically with 
the same number of steps per run. We chose the number 
of Monte Carlo steps (MCS) in each iteration (run) to 
increase linearly with the iteration number, r, specifically 
as r X 10^ MCS. This is a rather arbitrary choice which is 
necessary for comparison purposes. In fact, our proposal 
allows us to set the number of steps in the yet unexplored 
energy region as the criteria for moving on to the next 
run, which, in turn, allows better statistical control over 
the next estimate for the entropy. For a fair comparison 
between the two algorithms we choose /3o = for CFP 
and impose P{E) = in the histogram based method 
for E > iiQ. In this way the algorithms only explore the 
positive temperature region of the spectra. 

Figure 13 shows 7 histograms out of 35 runs of simu- 
lations on a 64 X 64 Ising model, with r x 10^ MCS per 
run, which was enough to reach the ground state with 
our algorithm. Several more runs are required using an 
histogram based method (lower panel of fig.|2l. It is im- 
portant to note that the full range of visited energies, 
~ (/io ~ Mr) I does not increase linearly with run number, 
r, because the spectral range added in each run shrinks 
as one approaches the ground state. In this respect the 
difference in performance of the two methods is quite re- 
markable. Histogram based methods do not provide an 
estimate for the entropy outside the previously visited 
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Figure 4: Full range of visited energies {Er — Eo)/N, as a 
function of the system size and number of iterations, r, for 
current method (left panel) and for histogram based recursion 
(right panel); in our moment based method r scales L = y/N 
whereas, with histograms, it scales with A*'. 



energy range. As a result, for energies in the previously 
unvisited range, H{E) decays as 

H{E) - exp [S{E)] - exp [-^riEr ~ E)] 

where Er is the lowest energy of the previous run and 
Pr = P{Er)- Our method uses uj{E) = firE which leaves 



H{E) - exp [S{E) - u{E)] - exp 



Since scales linearly with iV, the added energy range 
in each iteration scales differently in the two methods. 
In terms of energy per particle, the added energy range 
per run, Ae = AE/N, scales as 1/L = 7V~^/^ in the 
current proposal instead of {PrN) ^ in histogram based 
methods. In Figure^lwe plot the visited energy per par- 
ticle range after r iterations, for various system sizes, as 
a function of r/L for the CFP method and r/N for his- 
togram based method. The collapse of the curves for the 
various systems sizes shows that the number of runs re- 
quired to cover the same energy per particle range scales 
with with L — \/N in the current proposal (left panel) 
and TV in histogram based methods (right panel). 

The two methods perform similarly for small systems, 
but the advantage of CFP method becomes obvious for 
large N and no amount of fine tuning can disguise it. 

Nevertheless, other alternative schemes to Berg's re- 
cursion have already been proposed like Wang-Landau 
sampling |[2J| or the transition matrix method (2^. As 
will be seen shortly, the main advantage of our method is 
that, unlike these previous methods, it produces an an- 
alytic approximation to the microcanonical inverse tem- 
perature, (3{E), in an increasing energy range, right from 
the start of the simulation. That proves an asset in the 
implementation of procedures designed to overcome the 
slow down with system size that affects multicanonical 
simulations [lllTal2^. 




Figure 5: Tunneling time (scaled by (A'^ log 7V))^, for system 
sizes A'' = 48 X 48 and A'^ = 64 x 64, as an function of the 
recursion run divided by L which is proportional to the en- 
ergy per particle range, (a) CFP without optimization; (b) 
for an application of the optimization procedure only after an 
estimation of density of states in the chosen range of temper- 
atures; (c) our implementation of the optimization. The inset 
shows the histograms obtained for the sample with L — 64 
in several runs during the exploration of the energy spectrum 
with optimization. 



III. OPTIMIZATION OF TUNNELING TIMES 

In a recent publication, Trebst, Huse and Troyer 
showed that it is possible to decrease significantly the 
tunneling times of a multicanonical simulation, by aban- 
doning the requirement of a flat histogram |l8j |. Our 
procedure of construction of the statistical entropy is 
well suited to implement their optimization method, right 
from the start of the simulation, without having to con- 
struct an approximation to n{E) for the entire spectrum. 

The procedure proposed in |18| minimizes the average 
time required to span the gap between two fixed energies 
in the spectrum, and £'+. To achieve this purpose, 
one must distinguish each energy entry during the sim- 
ulation according to which of the two energies, E- or 
was visited last. We can thus measure separately 
n^{E), the number of visits to energy E occurring when 
the simulation has visited E- more recently than 
and n^(E) for the other way around. To minimize the 
tunneling times between the two energies i?_ and E^ 
one must choose an asymptotic energy distribution that 
satisfies 



H{E) cx 



dfjE) 
dE 



(7) 



where f{E 



n^(E)/H{E). The implementation pro- 
posed in [18j used the knowledge of the density of states, 
n{E) in the whole spectrum to measure the f{E), fol- 
lowing with a recursion procedure that converges to the 
asymptotic distribution which satisfies eq. [7| In our 
method of construction of the density of states there are, 
at each step, two energies which are the current bound- 
aries of the known density states, namely /Hq (that re- 
mains fixed) and /i^ (that changes with each run, r). 
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Therefore, by using E- = fir and E^ = /iq we can 
measure f{E). In the first run, for fii < Ei < /ip, we 
observe Hi{Ei) cx exp (—wi (£')), with loi{E) defined in 
eq. El rather than the optimal distribution H°p^{E) oc 
exp [In (d/i/dE')]. To converge to the optimal distri- 
bution (changing the weight alters f{E)), we use, in 
the next run, the geometric mean ^/HiH°p^ in the in- 
terval Hi < E < 01 : the weight factor becomes 
—uj2{E) + 4>2{E), where 



' 2 \dE 



for fii < E < ^0, 



with constant values of 02 (mi) for E < fii and 02 (mo) for 
E > HQ. Notice that the correction to the microcanoni- 
cal temperature, (3{E), is —dip2{E)/dE which is zero for 
E < Hi. This choice ensures the convergence to the cri- 
terion of eq. 13 in the range where the entropy is already 
known and where f{E) was measured, /ii < £^ < /iq, and 
gives a flat histogram in the region where the entropy 
was estimated by calculating moments of the Gaussian 
tail, i.e. ^-2 < E < fii- This procedure is iterated in the 
following runs with (j)r{E) defined as. 



cj>r{E) 



In 



d}r-l 

dE 



Mr-l < E < flQ 



with constant values (j)r{Hr-i) for E < fir-i and ^^(mo) 
for E > HQ. To extract the numerical derivative of f{E), 
avoiding the difficulties of the fluctuations in histogram 
entries, we use the natural scale afforded by cr^, and cal- 
culate df / E as 



dE 



{f{E) {E-Hp))p 



In figureElwe plot the average tunneling time, divided 
by (A^lnA^)^, for two different system sizes, as an func- 
tion of the run number of the multicanonical recursion 
divided by i = ^/N ; simulations with the same horizon- 
tal coordinate correspond to the same energy per particle 
range. Curves (a) correspond to the situation without op- 
timization, and the average tunneling times vary faster 
than (A^lniV)^; the two sizes do not collapse to a single 
curve. In case (c), with the optimization carried out while 
the density of states is being determined, the curves for 
the two system sizes track each other. If the optimiza- 
tion correction is only performed after full exploration 
of the energy spectrum, curve (b), there is no further 
gain in tunneling time, as the average tunneling times 
of (b) merge with curve (c). This observation clearly 
supports our suggestion that optimization can be imple- 
mented while the density of states is being constructed. 
In this fashion, it not only reduces tunneling times when 
n{E) is known, but als o sp eeds up the actual calculation 
of n{E). As found in [ig, the inset shows a strong sig- 
nature of the optimizing procedure in the critical energy 
where the diffusivity is low. 



IV. CLUSTER DYNAMICS 

An alternative way to improve the efficiency of mul- 
ticanonical simulations consists in changing the specific 
dynamics of the Markov chain, i.e., the algorithms used 
to propose and accept configuration moves. 

In canonical ensemble simulations, cluster update al- 
gorithms like Wolff's 0, Swendsen-Wang |^ or Nie- 
dermeyer's Q have proved very effective in overcoming 
critical slowing down of correlations. Several propos- 
als have been presented to generalize cluster update ap- 
proaches in multicanonical ensemble simulations, either 
using spin-bond representations of the partition function, 
,24, 25, 26, 27], or cluster building algorithms based on 
alternative ways of computing the microcanonical tem- 
perature, i3{E) mm. 

Wolff's cluster algorithm provides a clever way of 
growing a cluster of parallel spins which can be flipped 
with probability 1, and still maintain the required Boltz- 
mann asymptotic distribution. This remarkable possibil- 
ity is intimately related to the fact that lo{E) is linear in 
energy, uj[E) = f3E, in a Metropolis simulation. 

In a multicanonical simulation, each step of the corre- 
sponding Markov chain occurs with the same probability 
as that of a canonical ensemble simulation, with an effec- 
tive temperature Pi chosen as (3i — (3{Ei), Ei being the 
energy of the current configuration; hence the designation 
"Multicanonical" . With this in mind, the simplest way of 
implementing cluster dynamics in a multicanonical sim- 
ulation is to use Pi = l3{Ei) to grow a cluster exactly 
as proposed in Wolff's algorithm. However, since the re- 
verse path implies a different value of /3, (3j = (3{Ej), 
where Ej is the energy of the next configuration in the 
chain, we must include an acceptance probability to en- 
sure detailed balance. 

In figure El we illustrate a move involving the flipping 
of four spins (labelled 1 to 4) on the left, and the reverse 
move on the right. The site marked with the number 
1 has been chosen with uniform probability. If a bond 
connects spin 1 to a neighboring spin parallel to 1 it is 
added to the cluster with a probability pi and rejected 
with probability 1—pi. This step is then repeated for the 
neighbours of the initial spin which were added to the 
cluster, until the process stops and there are no further 
bonds that can be aggregated to the cluster. The proba- 
bility of generating a cluster with Ha accepted bonds, in 
which rir bonds to spins parallel to the initial one were 
inspected and rejected, is given by: 

G,^/=p^(i-p,r'^- 

It is important to note that includes a number of re- 
jected bonds that now link spins inside the cluster (like 
the bond from spin 1 to spin 4 in fig. jSJ. We write 
rir = Hp + Hf, where Hp counts the number of such bonds 
and n f is the number of bonds from spins in the cluster 
to parallel spins outside the cluster. This distinction is 
important when considering the reverse move. 



6 




• • • 



Figure 6: This scheme shows, on the left, how, starting from 
a given spin (marked as 1), the bonds to neighbouring spins 
are inspected and marked according to the algorithm. On the 
right, are shown the resulting configuration and the way the 
initial state can be reached from it. 




The difference in energy between the final and initial 
configuration is determined only by the frontier of the 
cluster; it is given as A_B — ~2J{nd — rif) where Ud is 
the number of bonds to spins opposite to the spins in 
the cluster. Let us now consider the reverse move, which 
requires us to select the same cluster in the same order 
with the spins now reversed we respect to the original 
state. 

Referring once again to figure one can see that the 
ria bonds that were accepted in the direct move (left) 
with probability Pi must be accepted in the reverse move 
(right), each with probability pj] the bonds to oppo- 
site spins in the direct move now connect to spins parallel 
to those in the cluster and must be rejected probability 
(1 — pj); the bonds to spins that were rejected in the 
direct move (n/) are, in the reverse process, bonds to 
opposite spins and are, hence, rejected with probability 
1; finally the Up bonds that were rejected but link spins 
inside the cluster, must now also be rejected. In other 
words, for the direct move 

G.^/=p^(i-P.r'(i-p.r^ (8) 

while, for the reverse process, 

Gf^,^p-^{i-p,r^{i-p,r^. (9) 

Wolff's algorithm corresponds to choosing pi = pj = 1 — 
exp(— 2 J/?) which implies that 

= {1 - pif^-^'f ^ e"^"" . 

The detailed balance condition for Boltzmann's equilib- 
rium distribution is obtained for an acceptance probabil- 
ity of 1 for flipping the cluster. To ensure an asymptotic 
distribution proportional to exp [—S {E)] , the detailed 
balance requires an acceptance probability given by 

A.^/=min(^l,g^^e-^^(^)). (10) 

If we choose pi — \ — e~'^^''' where f3i = (3{Ei), we 
find that this acceptance probability remains close to 1 
for most of the energy range (see fig. [TJ, falling only 



Figure 7: Acceptance ratio of the Wolff 's cluster algorithm, 
with a microcanonical /3{E), for the 2D ferromagnetic Ising 
model. In the paramagnetic phase the acceptance ratio grows 
to 1 with the increase of system size; in the ferromagnetic 
phase the acceptance ratio tends to a finite value, smaller 
than 1, as A'' grows. 

at very low temperatures. This behaviour is in strong 
contrast with the one for single spin fiip dynamics, where 
the acceptance rate is only 1 at the maximum of the 
density of states. If the inverse temperatures of the initial 
and final states are close, f3i « /3j, the acceptance rate 
becomes 

A-.f « min(l,(l-p,)"'~"''e-^^(^)) (11) 
= min(l,e^'^^e"^^(^)) (12) 

where we used AE ~ —2 J {rid — "■/) to obtain the second 
expression. For small energy differences, A5 « /3iAE 
and the acceptance rate becomes close to unity. 

In the case of the 2D ferromagnetic Ising model, the 
average excitation energy, AE, is of 0{Vn') for energies 
below p,c {p-c = Pp^ where /3c is the inverse critical tem- 
perature), of 0(1) for energies above p,c, with a crossover 
between these regimes in the neighborhood of (ic (fig-IHll- 
This behaviour reflects a huge difference with respect to 
single spin flip dynamics (SSF), where AE is always of 
0(1). A similar behaviour exists for the number of spins, 
n.u that are inspected in each call of the Markov chain: 
Uy ^ 0(1) for E > He and n„ ^ 0{N) for E < pc- This 
two regime behaviour introduces an additional complex- 
ity in this method. In particular, the tunneling time 
measured in Markov chains calls, no longer scales as the 
computational time with system size, since Markov chain 
calls can take a computational time of order O(A^). We 
therefore redefine the time scale so that a Markov chain 
call in a state of energy E corresponds to a time span 
of ny{E), the average number of inspected spins in the 
cluster buildup process. 

We now consider the system's coarse grained random 
walk in energy space. When the energy is close to E, the 
mean square energy change in M Markov chain calls is 
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Figure 8: Two regime behavior of the WolfT's cluster algo- 
rithm for the 2D ferromagnetic Ising model, with a micro- 
canonical temperature P{E). The left panel represents the 
equal energy average of the fraction of visited spins, n„, during 
the cluster growing. In the ferromagnetic phase the number 
of visited spins is of the order of 0{N) and in the paramag- 
netic phase is of 0(1). The right panel it is represents the 
equal energy average of the excitations, A_E. In the param- 
agnetic phase (AE^) becomes independent of the system 
size, while in the ferromagnetic phase it scales as \/iV. 



and this occurs in a time tm — M x ny{E). With time 
measured in this way, the diffusivity of this random walk 
in energy space is 



D{E) cx 



n,{E) ■ 



(13) 



On the other hand, the probability that the system is at 
energy E is proportional to 



H,{E)^H{E) X n,{E), 



(14) 



since each visit to energy E lasts a time ny{E). The 
probability current is given, quite generally, by 



3 = p{E)V{E) - — {D{E)p{E)) 



where, in equilibrium, j — 0, and 



p{E) - po{E) 



HyjE) 

J dEHy{E)' 



V{E) is a bias field, in general non-zero, which, together 
with D{E), determines the equilibrium distribution. 

It can be shown |0 [s^ that the tunneling times of 
a random walk with a given diffusivity, D{E), can be 
minimized with an optimal choice of bias field V{E). The 
corresponding equilibrium distribution is given by 



Po{E) 



1 



^D{E) 

Using eg nations 1 1 31 and 1 1 41 this condition becomes 

H{E) cx ^ = 

^{AE^)ny{E) 



(15) 




Figure 9: Results of the CFP recursion to construct an his- 
togram given by eq. 1151 (a) Histogram of the simulation for 
several steps of the CFP recursion; (b) Acceptance ratio; (c) 
{H{E) X y^njE){AE^) ; (d) Tunneling time (scaled by N'^) 
versus the run of the recursion (scaled by L). 



This variation of log {H{E)), relative to a flat histogram, 
is of order 0(log(iV)), and, therefore, histograms remain 
broad, covering the entire spectrum, but are no longer 
flat as shown in panel (a) of figure |51 On panel (d) it is 
shown that the tunneling time for this simulation scales 
has N'^ as expected from a simple diffusion. 

The optimization procedure presented in reference |18| , 
in the context of N-fold Way dynamics, is closely related 
to this one (but not identical) and leads to a choice 



H{E) cx 



df 



ny{E) dE 



where f{E) was defined above. We find only a marginal 
improvement in tunneling times with respect to the case 
of a histogram defined by eq. El These two procedures 
will be compared in another publication jsfj . 



V. CONCLUSION 

We have proposed a new method to build the density 
of states in an multicanonical simulation. The method is 
based on the calculation of moments of the energy dis- 
tribution. It avoids the use of histograms and can just 
as easily be implemented for continuous as for discrete 
systems. It leads to a piecewise analytic approximation 
to the microcanonical inverse temperature P{E). In any 
stage of the simulation there are two well defined ener- 
gies, E^ and E+, that limit the range in which f3{E) is 
known. Therefore the method can be applied without 
difficulty to a predefined temperature range such as a 
neighborhood of a critical temperature. 

We have also demonstrated the usefulness of this 
method in the implementation of various optimization 
schemes that render the simulation more efficient. In fig- 
ure ^1 we sum up the results we obtained for the scaling 
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Figure 10: Scaling of tunneling time with system size in vari- 
ous broad histogram methods: (SSF) straight forward multi- 
canonical simulation with single spin flip dynamics; ( SSFopt) 
optimized single spin flip dynamics; (clusteroi,) optimized bias 
for the measured D{E); (clusteropt) optimized ensemble with 
cluster dynamics; (cluster/;i) flat histogram with cluster dy- 
namics. 

of tunneling times with the system size. In general the 
scaling of the average tunneling time is r ^ ]\['^+^- . In 
a straightforward multicanonical simulation with a SSF 



dynamics, z = 0.39. Using the optimization procedure 
of j3| we confirm that r ^ (NlnN)"^. For a general- 
ization of Wolff's cluster method for the multicanonical 
ensemble we found a biased random walk in energy with 
z — 0.82. We proposed a new method for reducing tun- 
neling time of cluster update simulations which adjusts 
the bias of the random walk in energy space. In this 
case, (clusterob) and also for optimized ensemble simula- 
tion with cluster algorithm, proposed in (clusteropt), 
the results are compatible with z = 0. In the (clusteroh) 
method, however, one avoids the necessity to calculate 
of the derivative of f{E), required for the (clusteropt) 
method of [l^ . In terms of actual computer time, we 
also found that the amplitudes of the scaling laws are 
considerably smaller for the optimized cluster methods, 
than for SSF dynamics. 
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